En este trabajo realizaremos un intento de predecir las hipotecas de la provincia de Cáceres durante en el año 2019 empleando un modelo ARIMA en el contexto de las series temporales.
Importamos todas las librerías y asignamos el directorio de trabajo, además de cargar el fichero csv que contiene nuestros datos.
library(ggplot2)
library(plotly)
library(tseries)
library(MASS)
library(forecast)
library(lmtest)
library(caschrono)
library(lubridate)
library(timeDate)
library(tsoutliers)
library(dplyr)
# Se limpia el espacio de trabajo:
rm(list = ls())
# Se asigna el directorio donde están los datos
setwd("C:/Users/Javier/Documents/masterAFI/16. Analisis de series temporales/01. Series temporales/javier_herraez_albarran_series_temporales/")
# setwd("C:/Users/jherraez/Documents/masterAFI/16. Analisis de series temporales/01. Series temporales/practica/")
datos <- read.csv("_11_CACERES.csv")
Dividimos en train y en test, y convertimos los datos en un objeto de tipo time series para poder aplicar sobre ellos las funciones de R.
datos$FECHA = as.Date(datos$FECHA,format='%d/%m/%Y')
datos.train <- subset(datos, FECHA<=as.Date('01/12/2018',format='%d/%m/%Y'))
datos.test <- subset(datos, FECHA>as.Date('01/12/2018',format='%d/%m/%Y'))
datos.train.ts <- as.ts(datos.train$TARGET, frequency = 12)
datos.test.ts <- as.ts(datos.test$TARGET, frequency=12)
Elaboramos un gráfico que nos permita ver la evolución de las hipotecas en los últimos años.
names(datos.train)
## [1] "FECHA" "TARGET"
graficoInicial <- ggplot(aes(x= FECHA, y = TARGET), data = datos.train) +
geom_line(color = '#d84519', size = 1) +
xlab('FECHA') + ylab('Matriculaciones')
ggplotly(graficoInicial)
Vemos como se han reducido considerablemente las hipotecas en estos años más próximos, aunque se puede seguir intuyendo algo de periodicidad.
Evaluamos si es necesario transformar la serie para hacerla estacionaria en varianza. Para ello realizamos un boxcox.
box_cox <- boxcox(TARGET ~ FECHA,
data = datos.train,
lambda = c(0, 0.5, 1))
lambda <- box_cox$x[which.max(box_cox$y)]
lambda
## [1] 0.1616162
Como lambda toma un valor cercano a cero, tomamos logaritmos.
datos.train$log_target=log(datos.train$TARGET)
datos.test$log_target=log(datos.test$TARGET)
datos.train.ts <- as.ts(datos.train$log_target, frequency = 12)
datos.test.ts <- as.ts(datos.test$log_target, frequency=12)
Probamos a hacer el test de Dickey-Fuller aunque no confiamos demasiado en su resultado, ya que es poco exigente.
adf.test(datos.train.ts)
##
## Augmented Dickey-Fuller Test
##
## data: datos.train.ts
## Dickey-Fuller = -1.6659, Lag order = 5, p-value = 0.7159
## alternative hypothesis: stationary
Todo parece indicar que va a ser necesaria una diferencia regular, pero nosotros ajustamos directamente sin nada.
Observamos los gráficos f.a.s. y la f.a.p. para realizar el primer ajuste.
acf(datos.train.ts, lag.max = 48, xlab = "Retardo",
main= "Función de autocorrelación simple")
pacf(datos.train.ts, lag.max = 48, xlab = "Retardo",
main = "Función de autocorrelación parcial")
Lo habitual es empezar con una estructura AR(1). La fuerte correlación de orden 1 en la f.a.p parece reforzar esta decisión, por lo tanto, se propone un ARIMA(1,0,0). Vemos los coeficientes de las variables que obtenemos al realizar este primer modelo.
ajuste1 <- Arima(datos.train.ts,
order = c(1,0,0),
seasonal = list(order = c(0,0,0), period = 12),
method = "ML")
ajuste1
## Series: datos.train.ts
## ARIMA(1,0,0) with non-zero mean
##
## Coefficients:
## ar1 mean
## 0.9059 5.5001
## s.e. 0.0303 0.2199
##
## sigma^2 = 0.09146: log likelihood = -42.67
## AIC=91.34 AICc=91.47 BIC=101.11
coeftest(ajuste1)
##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 0.905855 0.030348 29.849 < 2.2e-16 ***
## intercept 5.500117 0.219941 25.007 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Probamos también a ver las correlaciones de estas.
cor.arma(ajuste1)
## ar1 intercept
## ar1 1.00000000 -0.02691084
## intercept -0.02691084 1.00000000
Vemos que en valor absoluto estarían por debajo de 0.8 por lo que no las consideremos correladas.
Mediante el test de LJung-Box estudiamos si hay ruido blanco.
Box.test.2(residuals(ajuste1),
nlag = c(6,12,18,24,30,36,42,48),
type="Ljung-Box")
## Retard p-value
## [1,] 6 6.49e-06
## [2,] 12 5.80e-07
## [3,] 18 2.54e-06
## [4,] 24 1.10e-07
## [5,] 30 6.80e-07
## [6,] 36 8.00e-08
## [7,] 42 4.30e-07
## [8,] 48 3.00e-07
Buscamos obtener unos p-valores mayores a 0.05, por lo que como vemos no tenemos ruido blanco.
Seguimos ajustando el modelo, así que observamos los gráficos f.a.s. y la f.a.p. para realizar un segundo ajuste.
acf(ajuste1$residuals, lag.max = 48, xlab = "Retardo",
main= "Función de autocorrelación simple")
pacf(ajuste1$residuals, lag.max = 48, xlab = "Retardo",
main = "Función de autocorrelación parcial")
Según lo que observamos en las gráficas, como destaca la línea del 1 en la gráfica de autorrelación simple, ajustamos un modelo ARIMA(1,0,1).
ajuste2 <- Arima(datos.train.ts,
order = c(1,0,1),
seasonal = list(order = c(0,0,0), period = 12),
method = "ML")
ajuste2
## Series: datos.train.ts
## ARIMA(1,0,1) with non-zero mean
##
## Coefficients:
## ar1 ma1 mean
## 0.9947 -0.7254 5.5238
## s.e. 0.0054 0.0434 0.5351
##
## sigma^2 = 0.06019: log likelihood = -2.52
## AIC=13.04 AICc=13.26 BIC=26.07
coeftest(ajuste2)
##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 0.9947456 0.0054474 182.610 < 2.2e-16 ***
## ma1 -0.7254036 0.0434252 -16.705 < 2.2e-16 ***
## intercept 5.5238319 0.5351063 10.323 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Vemos que el coeficiente del AR1 + 2 * Std. Error pasaría del valor 1, así que es claro que necesitamos diferenciarlo. Pasamos entonces a un modelo ARIMA(0,1,1).
ajuste3 <- Arima(datos.train.ts,
order = c(0,1,1),
seasonal = list(order = c(0,0,0), period = 12),
method = "ML")
ajuste3
## Series: datos.train.ts
## ARIMA(0,1,1)
##
## Coefficients:
## ma1
## -0.7276
## s.e. 0.0426
##
## sigma^2 = 0.05974: log likelihood = -1.79
## AIC=7.59 AICc=7.65 BIC=14.09
coeftest(ajuste3)
##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ma1 -0.727557 0.042638 -17.064 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Hacemos el test de Ljung-Box
Box.test.2(residuals(ajuste3),
nlag = c(6,12,18,24,30,36,42,48),
type="Ljung-Box")
## Retard p-value
## [1,] 6 0.67336896
## [2,] 12 0.23278469
## [3,] 18 0.25043194
## [4,] 24 0.04279890
## [5,] 30 0.09096946
## [6,] 36 0.03682012
## [7,] 42 0.04956558
## [8,] 48 0.07392056
Vemos que ya casi todos los valores están por encima del 0.05, aún así seguimos ajustando.
acf(ajuste3$residuals, lag.max = 48, xlab = "Retardo",
main= "Función de autocorrelación simple")
pacf(ajuste3$residuals, lag.max = 48, xlab = "Retardo",
main = "Función de autocorrelación parcial")
Según las gráficas, vemos que podemos proponer un SARIMA(0,1,1)x(1,0,0)[12]. Construyamos este modelo..
ajuste4 <- Arima(datos.train.ts,
order = c(0,1,1),
seasonal = list(order = c(1,0,0), period = 12),
method = "ML")
ajuste4
## Series: datos.train.ts
## ARIMA(0,1,1)(1,0,0)[12]
##
## Coefficients:
## ma1 sar1
## -0.7430 0.2382
## s.e. 0.0445 0.0729
##
## sigma^2 = 0.05668: log likelihood = 3.35
## AIC=-0.71 AICc=-0.58 BIC=9.05
coeftest(ajuste4)
##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ma1 -0.743005 0.044523 -16.6881 < 2.2e-16 ***
## sar1 0.238164 0.072912 3.2664 0.001089 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Medimos también correlaciones.
cor.arma(ajuste4)
## ma1 sar1
## ma1 1.0000000 -0.1216682
## sar1 -0.1216682 1.0000000
Y volvemos a hacer el test de Ljung-Box para ver si ya obtenemos ruido blanco.
Box.test.2(residuals(ajuste4),
nlag = c(6,12,18,24,30,36,42,48),
type="Ljung-Box")
## Retard p-value
## [1,] 6 0.7123801
## [2,] 12 0.9055981
## [3,] 18 0.8138609
## [4,] 24 0.5663195
## [5,] 30 0.7499105
## [6,] 36 0.6405039
## [7,] 42 0.6366659
## [8,] 48 0.7534805
Vemos que tenemos valores por encima de 0.05 para todos, así que parece que podemos decir que tenemos ruido blanco.
Crearemos una función para crear variables explicativas en series mensuales. Esta función vista en clase añade alguna variante para añadir festivos propios de Cáceres.
calculoExplicativasCalendarioCaceres <- function(variableFecha, domingoYFestivosJuntos){
#######################################
# Creación de todas las fechas #
#######################################
if (month(max(variableFecha)) %in% c(1,3,5,7,8,10,12)) {
diasHastaFinMes <- 30
} else if (month(max(variableFecha)) %in% c(4,6,9,11)) {
diasHastaFinMes <- 29
} else if (year(max(variableFecha))%%4==0) {
diasHastaFinMes <- 28
} else {diasHastaFinMes <- 27}
todasLasFechas <- data.frame(fechas=seq(min(variableFecha),
max(variableFecha)+diasHastaFinMes,
by="days"))
#######################################
# Cálculo de la Semana Santa #
#######################################
domingoResurrecion <- as.Date(Easter(year(min(variableFecha)):year(max(variableFecha))))
lunesPascua <- domingoResurrecion + 1
sabadoSanto <- domingoResurrecion - 1
viernesSanto <- domingoResurrecion - 2
juevesSanto <- domingoResurrecion - 3
# Se unen y ordenan todos los días que forman la Semana Santa
semanaSanta <- sort(c(juevesSanto, viernesSanto, sabadoSanto, domingoResurrecion, lunesPascua))
# Se pone en formato data.frame y se añade un indicador
semanaSanta <- data.frame(fechas=semanaSanta, semanaSanta=rep(1,length(semanaSanta)))
# Se añaden a la tabla maestra de fechas
todasLasFechas_2 <- merge(x = todasLasFechas, y = semanaSanta, by = "fechas", all.x = TRUE)
# Se reemplazan los NAs por 0
todasLasFechas_2$semanaSanta[is.na(todasLasFechas_2$semanaSanta)] <- 0
######################################
# Cálculo de la variable dt #
######################################
# 1. Definición de festivos:
############################
calendario <- todasLasFechas
calendario$diaSemana <- as.factor(wday(calendario$fecha))
calendario$diaMes <- as.factor(day(calendario$fecha))
calendario$mes <- as.factor(month(calendario$fecha))
calendario$anyo <- as.factor(year(calendario$fecha))
calendario$p_01ene <- ifelse(calendario$diaMes==1 & calendario$mes==1, 1, 0)
calendario$p_06ene <- ifelse(calendario$diaMes==6 & calendario$mes==1, 1, 0)
calendario$p_19mar <- ifelse(calendario$diaMes==19 & calendario$mes==3, 1, 0)
calendario$p_01may <- ifelse(calendario$diaMes==1 & calendario$mes==5, 1, 0)
calendario$p_15ago <- ifelse(calendario$diaMes==15 & calendario$mes==8, 1, 0)
calendario$p_12oct <- ifelse(calendario$diaMes==12 & calendario$mes==10,1, 0)
calendario$p_01nov <- ifelse(calendario$diaMes==1 & calendario$mes==11, 1 ,0)
calendario$p_06dic <- ifelse(calendario$diaMes==6 & calendario$mes==12, 1 ,0)
calendario$p_08dic <- ifelse(calendario$diaMes==8 & calendario$mes==12, 1 ,0)
calendario$p_25dic <- ifelse(calendario$diaMes==25 & calendario$mes==12, 1 ,0)
# Fiestas regionales: día de Extremadura
calendario$p_08sep <- ifelse(calendario$diaMes==8 & calendario$mes==9, 1 ,0)
# Fiestas locales: San Jorge
calendario$p_23abr <- ifelse(calendario$diaMes==23 & calendario$mes==4, 1 ,0)
calendario$festivo <- rowSums(subset(calendario, select=p_01ene:p_25dic))
# La definición de la variable dt varía según la opción domingoYFestivosJuntos.
if (domingoYFestivosJuntos==0){
calendario$sabado <- ifelse(calendario$diaSemana==7, 1 ,0)
calendario$domingo <- ifelse(calendario$diaSemana==1, 1 ,0)
# Días laborables: todos menos sábados y domingos
calendario$laborable <- 1-calendario$sabado-calendario$domingo
} else {
calendario$sabado <- ifelse(calendario$diaSemana==7, 1 ,0)
# Domingo=1 si domingo=1 o festivo=1
calendario$domingo <- ifelse(calendario$domingo==1 | calendario$festivo==1, 1 ,0)
# Días laborables: todos menos sábados y domingos/festivos
calendario$laborable <- 1-calendario$sabado-calendario$domingo
}
# 2. Definición de variable dt:
###############################
# Se filtran las columnas de inter?s y se añade la Semana Santa
calendario_2 <- calendario[, c("fechas", "mes", "anyo", "sabado", "domingo", "laborable", "festivo")]
todasLasFechasFinal <- merge(x = todasLasFechas_2, y = calendario_2,
by = "fechas", all.x = TRUE)
# Agregamos la serie a nivel a?o-mes
calendarioAnyoMes <- aggregate(todasLasFechasFinal[,c("sabado","domingo",
"laborable", "semanaSanta", "festivo")],
by=list(mes=todasLasFechasFinal$mes,
anyo=todasLasFechasFinal$anyo),
"sum")
# Se calcula la variable dt:
calendarioAnyoMes$dt <- calendarioAnyoMes$laborable-(5/2)*(calendarioAnyoMes$sabado+calendarioAnyoMes$domingo)
######################################
# Cálculo de años bisiestos #
######################################
calendarioAnyoMes$anyoNum <- as.numeric(levels(calendarioAnyoMes$anyo))[calendarioAnyoMes$anyo]
calendarioAnyoMes$bisiesto <- ifelse(calendarioAnyoMes$mes==2 &(calendarioAnyoMes$anyoNum %% 4)==0, 1 ,0)
#######################################################
# Tabla final con explicativas de calendario #
#######################################################
if (domingoYFestivosJuntos==0){
explicativasCalendario <- cbind(fecha=variableFecha, calendarioAnyoMes[, c("semanaSanta", "dt", "bisiesto", "festivo")])
} else {
explicativasCalendario <- cbind(fecha=variableFecha, calendarioAnyoMes[, c("semanaSanta", "dt", "bisiesto")])
}
return(explicativasCalendario)
}
Llamamos a la función que hemos creado con la opción de que domingos y días festivos se cuenten por separado.
explicativasCalendarioTrain <- calculoExplicativasCalendarioCaceres(datos.train$FECHA, domingoYFestivosJuntos=0)
calendarioTrain <-
as.matrix(
explicativasCalendarioTrain[,c("semanaSanta", "dt", "bisiesto", "festivo")]
)
Añadimos al ajuste las nuevas variables explicativas.
ajuste4conCalendario <- Arima(datos.train.ts,
order = c(0,1,1),
seasonal = list(order = c(1,0,0), period = 12),
xreg = calendarioTrain,
method = "ML")
ajuste4conCalendario
## Series: datos.train.ts
## Regression with ARIMA(0,1,1)(1,0,0)[12] errors
##
## Coefficients:
## ma1 sar1 semanaSanta dt bisiesto festivo
## -0.7200 0.1969 -0.0116 0.0150 -0.0227 -0.0614
## s.e. 0.0483 0.0734 0.0118 0.0056 0.1049 0.0212
##
## sigma^2 = 0.0536: log likelihood = 10.88
## AIC=-7.76 AICc=-7.15 BIC=15.01
coeftest(ajuste4conCalendario)
##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ma1 -0.7199518 0.0483371 -14.8944 < 2.2e-16 ***
## sar1 0.1968911 0.0733748 2.6834 0.007289 **
## semanaSanta -0.0115610 0.0118108 -0.9788 0.327654
## dt 0.0150418 0.0055963 2.6878 0.007192 **
## bisiesto -0.0227074 0.1049431 -0.2164 0.828693
## festivo -0.0614426 0.0211852 -2.9003 0.003729 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Vemos que el año bisiesto tiene un p-valor bastante grande.Además, el valor que puede aportarnos tener un día más en febrero ya está contemplado con la variable dt. Por otro lado, podemos mantener la Semana Santa ya que le podemos ver sentido de negocio.
calendarioTrain <-
as.matrix(
explicativasCalendarioTrain[,c("semanaSanta", "dt", "festivo")]
)
ajuste5conCalendario <- Arima(datos.train.ts,
order = c(0,1,1),
seasonal = list(order = c(1,0,0), period = 12),
xreg = calendarioTrain,
method = "ML")
ajuste5conCalendario
## Series: datos.train.ts
## Regression with ARIMA(0,1,1)(1,0,0)[12] errors
##
## Coefficients:
## ma1 sar1 semanaSanta dt festivo
## -0.7203 0.1976 -0.0114 0.0150 -0.0608
## s.e. 0.0482 0.0733 0.0118 0.0056 0.0210
##
## sigma^2 = 0.05333: log likelihood = 10.86
## AIC=-9.71 AICc=-9.26 BIC=9.8
coeftest(ajuste5conCalendario)
##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ma1 -0.7202627 0.0482284 -14.9344 < 2.2e-16 ***
## sar1 0.1975910 0.0732627 2.6970 0.006996 **
## semanaSanta -0.0113950 0.0117906 -0.9664 0.333821
## dt 0.0150369 0.0055984 2.6859 0.007233 **
## festivo -0.0607887 0.0209877 -2.8964 0.003775 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Vemos correlaciones y vemos que ninguna variable parece tener correlación:
cor.arma(ajuste5conCalendario)
## ma1 sar1 semanaSanta dt festivo
## ma1 1.000000000 -0.1395694285 -0.02512203 -0.0060496156 -0.07048090
## sar1 -0.139569429 1.0000000000 -0.07245826 -0.0004302792 -0.03481278
## semanaSanta -0.025122031 -0.0724582600 1.00000000 0.0491736995 0.16027535
## dt -0.006049616 -0.0004302792 0.04917370 1.0000000000 0.00428523
## festivo -0.070480901 -0.0348127821 0.16027535 0.0042852305 1.00000000
Y el text de Ljung-Box para ruido blanco.
Box.test.2(residuals(ajuste5conCalendario),
nlag = c(6,12,18,24,30,36,42,48),
type="Ljung-Box")
## Retard p-value
## [1,] 6 0.7984990
## [2,] 12 0.9116931
## [3,] 18 0.8911888
## [4,] 24 0.7504270
## [5,] 30 0.9289244
## [6,] 36 0.8699718
## [7,] 42 0.8579565
## [8,] 48 0.9146346
Como resultado obtenemos que tenemos ruido blanco.
Pasamos a identificar outliers con la función locate.ouliers.
Identificamos outliers de los tipos:
listaOutliersTrain <- locate.outliers(ajuste5conCalendario$residuals,
pars = coefs2poly(ajuste5conCalendario),
types = c("AO", "LS", "TC"),
cval=3)
listaOutliersTrain$abststat=abs(listaOutliersTrain$tstat)
# Cruzamos con la tabla original para obtener la fecha
datos.train$ind <- as.numeric(rownames(datos.train))
listaOutliersTrainFecha <- merge(listaOutliersTrain, datos.train[,c("ind", "FECHA")], by = "ind")
arrange(listaOutliersTrainFecha,desc(listaOutliersTrainFecha$abststat))
## ind type coefhat tstat abststat FECHA
## 1 99 LS -0.4456849 -3.135621 3.135621 2011-03-01
## 2 126 TC -0.5369467 -3.119897 3.119897 2013-06-01
## 3 88 TC 0.5255352 3.053591 3.053591 2010-04-01
Nos han salido 2 outliers de tipo TC y 1 de tipo LS. Estos outliers los añadimos a las variables explicativas que ya teníamos.
outliersTrain <- outliers(c( "TC","LS","TC"), c(88, 99, 126))
outliersVariablesTrain <- outliers.effects(outliersTrain, length(ajuste5conCalendario$residuals))
calendarioMasOutliersTrain <- as.matrix(cbind(calendarioTrain, outliersVariablesTrain))
Ajustamos un nuevo modelo con todo ello.
ajuste5conCalendarioYOutliers <- Arima(datos.train.ts,
order = c(0,1,1),
seasonal = list(order = c(1,0,0), period = 12),
xreg = calendarioMasOutliersTrain,
method = "ML")
coeftest(ajuste5conCalendarioYOutliers)
##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ma1 -0.7938138 0.0404222 -19.6381 < 2.2e-16 ***
## sar1 0.1901647 0.0759103 2.5051 0.0122409 *
## semanaSanta -0.0141030 0.0114148 -1.2355 0.2166407
## dt 0.0139374 0.0054566 2.5542 0.0106425 *
## festivo -0.0633449 0.0197687 -3.2043 0.0013539 **
## TC88 0.4022273 0.1792464 2.2440 0.0248330 *
## LS99 -0.3986618 0.1380351 -2.8881 0.0038755 **
## TC126 -0.5954287 0.1754221 -3.3943 0.0006881 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Si miramos las correlaciones vemos que no existen valores altos:
cor.arma(ajuste5conCalendarioYOutliers)
## ma1 sar1 semanaSanta dt festivo
## ma1 1.000000000 -0.102755031 -0.02714986 0.005189159 -0.01367233
## sar1 -0.102755031 1.000000000 -0.03635183 -0.004683842 -0.01923234
## semanaSanta -0.027149858 -0.036351834 1.00000000 0.053407042 0.15537449
## dt 0.005189159 -0.004683842 0.05340704 1.000000000 0.00851406
## festivo -0.013672329 -0.019232340 0.15537449 0.008514060 1.00000000
## TC88 0.139030350 -0.008109016 -0.12474959 -0.013125127 0.05765664
## LS99 0.149116326 0.182599302 -0.07731113 -0.025581756 0.04218715
## TC126 0.164535310 0.134369279 0.03476912 0.074418665 0.06490375
## TC88 LS99 TC126
## ma1 0.139030350 0.14911633 0.16453531
## sar1 -0.008109016 0.18259930 0.13436928
## semanaSanta -0.124749593 -0.07731113 0.03476912
## dt -0.013125127 -0.02558176 0.07441866
## festivo 0.057656643 0.04218715 0.06490375
## TC88 1.000000000 0.27256859 0.02299857
## LS99 0.272568590 1.00000000 0.05621866
## TC126 0.022998565 0.05621866 1.00000000
Volvemos a realizar el test de Ljung-Box y vemos que nos vuelve a dar ruido blanco.
Box.test.2(residuals(ajuste5conCalendarioYOutliers),
nlag = c(6,12,18,24,30,36,42,48),
type="Ljung-Box")
## Retard p-value
## [1,] 6 0.8503906
## [2,] 12 0.8326711
## [3,] 18 0.8147331
## [4,] 24 0.7642923
## [5,] 30 0.9193775
## [6,] 36 0.9386808
## [7,] 42 0.9294358
## [8,] 48 0.9033291
Una vez ya obtenido nuestro modelo, vamos a pasar a intentar predecir el número de hipotecas que se concedieron en el año 2019.
Lo primero será añadir las variables explicativas del calendario a test.
explicativasCalendarioTest <- calculoExplicativasCalendarioCaceres(datos.test$FECHA,domingoYFestivosJuntos=0)
calendarioTest <- as.matrix(explicativasCalendarioTest[,c("semanaSanta", "dt", "festivo")])
Para los outliers tomaremos los del último año y los calcaremos en el año que vamos a predecir. Tras esto, juntamos las variables del calendario y las de los outliers.
outliersVariablesTest <- outliersVariablesTrain[(nrow(outliersVariablesTrain) - 11):nrow(outliersVariablesTrain),]
calendarioMasOutliersTest <- as.matrix(cbind(calendarioTest, outliersVariablesTest))
calendarioMasOutliersTest
## semanaSanta dt festivo TC88 LS99 TC126
## [1,] 0 3.0 2 3.927514e-15 1 3.022680e-09
## [2,] 0 0.0 0 2.749260e-15 1 2.115876e-09
## [3,] 0 -4.0 1 1.924482e-15 1 1.481113e-09
## [4,] 5 2.0 0 1.347137e-15 1 1.036779e-09
## [5,] 0 3.0 1 9.429961e-16 1 7.257455e-10
## [6,] 0 -5.0 0 6.600972e-16 1 5.080219e-10
## [7,] 0 3.0 0 4.620681e-16 1 3.556153e-10
## [8,] 0 -0.5 1 3.234477e-16 1 2.489307e-10
## [9,] 0 -1.5 0 2.264134e-16 1 1.742515e-10
## [10,] 0 3.0 1 1.584893e-16 1 1.219760e-10
## [11,] 0 -1.5 1 1.109425e-16 1 8.538323e-11
## [12,] 0 -0.5 3 7.765978e-17 1 5.976826e-11
Procedemos por lo tanto, a hacer la predicción de este nuevo año. Además, marcaremos los límites de confianza al 95%.
prediccion <- as.data.frame(predict(ajuste5conCalendarioYOutliers,
n.ahead=12,
newxreg=calendarioMasOutliersTest))
#Límites de confianza
U <- exp(prediccion$pred + 2*prediccion$se)
L <- exp(prediccion$pred - 2*prediccion$se)
datos.pred <- data.frame(FECHA = datos.test$FECHA,
Prediccion = exp(prediccion$pred + 0.5 * prediccion$se^2),
LimSup = U,
LimInf = L)
Para ver nuestro resultado, elaboramos una gráfica mostrando la predicción, los límites de confianza y el valor real de la serie para el año 2019.
datos.real.pred <- merge(datos, datos.pred, by = "FECHA", all.x = T)
datos.real.pred$REAL <- datos.real.pred$TARGET
datos.real.pred$TARGET[datos.real.pred$FECHA>as.Date('01/12/2018',format='%d/%m/%Y')] <- NA
grafico <- ggplot(data = datos.real.pred) +
geom_line(aes(x= FECHA, y = REAL), color = 'blue',
alpha = 0.4, size = 0.8) +
geom_line(aes(x= FECHA, y = TARGET), color = 'blue',
alpha = 0.8, size = 0.8) +
geom_line(aes(x= FECHA, y = Prediccion), color = 'darkred',
size = 1) +
geom_line(aes(x= FECHA, y = LimSup), color = 'orange',
size = 1) +
geom_line(aes(x= FECHA, y = LimInf), color = 'orange',
size = 1) +
xlab('FECHA') + ylab('Matriculaciones')
ggplotly(grafico)
Podemos, además, comparar nuestro modelo con el ajuste que se obtiene al hacer un auto arima. Veamos el resultado.
ajusteAutomatico <- auto.arima(datos.train.ts,
max.d=1, max.D=1,
max.p=2, max.P=2,
max.q=2, max.Q=2,
seasonal=TRUE,
ic="aic",
allowdrift=FALSE,
xreg=calendarioMasOutliersTrain,
stepwise=TRUE)
ajusteAutomatico
## Series: datos.train.ts
## Regression with ARIMA(2,0,0) errors
##
## Coefficients:
## ar1 ar2 intercept semanaSanta dt festivo TC88 LS99
## 0.3359 0.2129 6.1811 -0.0118 0.0144 -0.0766 0.3343 -1.2018
## s.e. 0.0730 0.0744 0.0558 0.0116 0.0051 0.0186 0.2375 0.0765
## TC126
## -0.6586
## s.e. 0.2284
##
## sigma^2 = 0.05661: log likelihood = 7.7
## AIC=4.6 AICc=5.82 BIC=37.18
El ajuste automático resuelve esta serie temporal con un ARIMA(2,0,0).
Vamos a calcular sus predicciones y a compararlas con las nuestras.
prediccion.auto <- as.data.frame(predict(ajusteAutomatico,
n.ahead=12,
newxreg=calendarioMasOutliersTest))
datos.pred.auto <- data.frame(FECHA = datos.test$FECHA,
Prediccion.auto = exp(prediccion.auto$pred + 0.5 * prediccion.auto$se^2))
datos.real.pred.auto <- merge(datos.real.pred, datos.pred.auto, by = "FECHA", all.x = T)
grafico <- ggplot(data = datos.real.pred.auto) +
geom_line(aes(x= FECHA, y = REAL), color = 'blue',
alpha = 0.4, size = 0.8) +
geom_line(aes(x= FECHA, y = TARGET), color = 'blue',
alpha = 0.8, size = 0.8) +
geom_line(aes(x= FECHA, y = Prediccion), color = 'darkred',
size = 1) +
geom_line(aes(x= FECHA, y = Prediccion.auto), color = 'orange',
size = 1) +
xlab('FECHA') + ylab('Matriculaciones') + ggtitle('Ajuste Manual vs Ajuste Automático')
ggplotly(grafico)
Vemos que nos arrojan unos resultados parecidos.